package fem2.material;

import inf.jlinalg.SolveFailedException;
import inf.text.ArrayFormat;

import org.junit.Test;

import fem2.Constraint;
import fem2.Edge;
import fem2.Element;
import fem2.Face;
import fem2.MaterialModel;
import fem2.Mesh;
import fem2.MeshUtilities;
import fem2.Model;
import fem2.Node;
import fem2.analysis.Analysis;
import fem2.analysis.NonlinearStaticAnalysis;
import fem2.element.BBarStructuralElement;
import fem2.element.EASStructuralElement;
import fem2.element.StructuralElement;
import fem2.enu.EASModeType;
import fem2.enu.EchoLevelType;
import fem2.enu.SolverType;
import fem2.enu.State;
import fem2.material.ExponentialSofteningRule;
import fem2.material.FluencyCriteria;
import fem2.material.IsotropicDamageMM;
import fem2.observer.StressObserver;

public class PatchTestIsotropicDamage {

	public MaterialModel createMaterial() {
		double E = 2;
		double nu = 0.3;
		double t = 1.0;
		State ss = State.PLANE_STRAIN;
		double ft = 0.91;
		// double Gf = 0.1;
		double e0 = ft / Math.sqrt(E);
		// double ef = Gf / ft;
		double ef = 100 * e0;
		FluencyCriteria softeningRule = new ExponentialSofteningRule(e0, ef);
		// return new IsotropicDamageMM(E, nu, 0, t, ss, damageLaw,
		// EquivalentStrainType.STRAIN_NORM);
		return new IsotropicDamageMM(E, nu, 0, t, ss, softeningRule, e0);
	}

	// public Mesh createMesh() { // the damage model works correctly on this
	// // mesh
	// Mesh mesh = new Mesh();
	//
	// Node n1 = mesh.addNode(new Node(0, 0));
	// Node n2 = mesh.addNode(new Node(1, 0));
	// Node n3 = mesh.addNode(new Node(1, 1));
	// Node n4 = mesh.addNode(new Node(0, 1));
	//
	// mesh.addFace(new Face(n1, n2, n3));
	// mesh.addFace(new Face(n1, n3, n4));
	//
	// mesh.setNumberOfFields(2);
	//
	// n1.setConstraint(new Constraint(false, false));
	// n4.setConstraint(new Constraint(false, true));
	// n2.setConstraint(new Constraint(true, false));
	//
	// return mesh;
	// }

	public Mesh createMesh() {
		/*
		 * mesh
		 */
		Mesh mesh = new Mesh();

		double l = 1;
		double h = 1;

		// double rx = Math.random();
		// double ry = Math.random();

		// double rx = 0.5;
		// double ry = 0.5;

		double rx = 0.17013914381185946;
		double ry = 0.5965964697460487;

		// Debugger.watch("rx = " + rx);
		// Debugger.watch("ry = " + ry);

		Node n1 = mesh.addNode(new Node(0, 0));
		Node n2 = mesh.addNode(new Node(l, 0));
		Node n3 = mesh.addNode(new Node(l, h));
		Node n4 = mesh.addNode(new Node(0, h));
		Node n5 = mesh.addNode(new Node(rx * l, ry * h));

		mesh.addFace(new Face(n1, n2, n5));
		mesh.addFace(new Face(n2, n3, n5));
		mesh.addFace(new Face(n3, n4, n5));
		mesh.addFace(new Face(n4, n1, n5));

		mesh.addEdge(new Edge(n2, n3));

		mesh.setNumberOfFields(2);

		/*
		 * constraint
		 */
		Constraint c1 = new Constraint(false, false);
		Constraint c2 = new Constraint(true, false);
		Constraint c3 = new Constraint(false, true);
		n1.setConstraint(c1);
		n2.setConstraint(c2);
		n4.setConstraint(c3);

		return mesh;
	}

	// public MaterialModel createMaterial() {
	// double E = 10;
	// double nu = 0.2;
	// double t = 1.0;
	// State ss = State.PLANE_STRAIN;
	// double e0 = 0.5 * Math.sqrt(E);
	// double ef = 1.2 * Math.sqrt(E);
	// FluencyCriteria softeningRule = new ExponentialSofteningRule(e0, ef);
	// return new IsotropicDamageMM(E, nu, 0, t, ss, softeningRule);
	// }

	public Model createModel(Mesh mesh, int elementType) {
		Model M = new Model(mesh);
		MaterialModel mm = createMaterial();

		for (int i = 0; i < mesh.countFaces(); i++) {
			Element e;

			switch (elementType) {
			case 0:
				e = new StructuralElement(mesh.getFace(i), mm);
				break;
			case 1:
				e = new BBarStructuralElement(mesh.getFace(i), mm);
				break;
			case 2:
				e = new EASStructuralElement(mesh.getFace(i), mm, EASModeType.Q1E4);
				break;
			default:
				throw new Error("invalid element type");
			}

			M.addElement(e);

			e.setEchoLevel(EchoLevelType.DEBUG, EchoLevelType.AUX2, EchoLevelType.OUTPUT);
		}

		return M;
	}

	@Test
	public void testSolve1() throws SolveFailedException {
		// Mesh mesh = new PatchTestT3(0).createMesh(); // does not work stably,
		// need check

		Mesh mesh = createMesh();
		Model model = createModel(mesh, 0);

		double du = 2;
		Constraint c1 = new Constraint(false, false);
		c1.setValue(0, du);
		Constraint c2 = new Constraint(false, true);
		c2.setValue(0, du);

		mesh.setConstraint(c2, MeshUtilities.seekNodesOnSurface(mesh, 1, 0, -1));
		MeshUtilities.seekNode(mesh, 1, 0).setConstraint(c1);

		// Analysis an = new LinearStaticAnalysis(model);
		Analysis an = new NonlinearStaticAnalysis(model, 10);
		an.getStrategy().setSolverType(SolverType.BAITSCH_SOLVER_GENERIC);
		an.getStrategy().addObserver(new StressObserver(model));
		an.getStrategy().setEchoLevel(EchoLevelType.OUTPUT, EchoLevelType.DEBUG);
		an.run();

		// double[][] Uh = new double[][] { { 0, 0 }, { 2.0, 0 }, { 2.0,
		// -0.967228028754730 },
		// { 0, -0.967228028754730 } };
		for (int i = 0; i < 4; i++) {
			double[] uHatLocal = model.getMesh().getNode(i).getUHatLocal();
			System.out.println(ArrayFormat.format(uHatLocal));
		}
		// for (int i = 0; i < 4; i++) {
		// double[] uHatLocal = model.getMesh().getNode(i).getUHatLocal();
		// AssertJUnit.assertEquals(Uh[i][0], uHatLocal[0], 1e-10);
		// AssertJUnit.assertEquals(Uh[i][1], uHatLocal[1], 1e-10);
		// }

		// IVector epsilon = ((StructuralElement)
		// m.getElement(0)).computeTotalStrain(0);
		// IVector sigma = ((StructuralElement) m.getElement(0)).getStress(0);
		// System.out.println("strain = " + epsilon);
		// System.out.println("stress = " + sigma);
	}
}
